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This work deals in evaluating the performance 
of various filters reported in literature and their 
corresponding higher order filters developed at IIT- 
Kanpur, Initially, a sensitivity analysis had been 
carried out for various filter pa ram eters using fan- 
beam convolution-backprojecti on algorithm for various 
simulated bubbly and annular flow data. The projection 
data for air— water bubbly flov/ is also processed to 
obtain the density distribution. Finally, the parallel- 
beam con volution-backprojecti on algorithm has been tested 
for various simulated bubbly and annular flow data, 

A comparative study of the relative performance of the 
linear filters (and their corresponding higher order 
filters), has been carried out with the simulated data. 
The results indicate an increase in both accuracy and 
smoothness if the higher order filtersi'‘are used. 



CHAPTER _ 
I MFRO DUCTION 


The algorithms used in present day CAT scanners 
are generally based on the convolution-backprojection 
scheme suggested by Ramachandran and Lakshminarayanan [1], 
This algorithm requires measured data in the form of line 
integrals of the spatial distribution of the physical 
property of interest Convolution~backprojection is so 
named because of the two basic steps involved in this 
algorithm: (a) measured data (called projection data), 
which has been collected by the scanner, is convolved 
with a filter function, and (b) the filtered data is 
backpro jected to form the image. The choice of these 
filters depends upon the user, A class of linear filters 
(in frequency domain) for the parallel beam geometry has 
been suggested by Lewitt [2] , The accuracy of recons- 
truction depends on the judicious choice of filter appro- 
priate for the function being investigated. This selection 
process (of appropriate filter) puts a severe limitation 
on the methodology which claims independence of a priori 
information re^.ding the nature cf function being investi- 
gated. 

The present study is an effort towaicls evaluating 
the performance of some existing filters, and their corres- 
ponding higher order filters developed at IIT Kanpur, The 



algorithm was initially tested with simulated data repre- 
senting bubbly and annular flovjs , The air-water flow 
data [3,4,5} was processed with linear filter to get 
density maps for the study of Ref. [^'j , Reconstruction 
of the simulateddata shows encouraging results for the 
higher order filters. 



CHAPTa, ^2 
mE U/ffNARlES 

In this chapter, we discuss the computational and 
mathematical procedures underlying the data collection, 
image reconstruction, and irrage display in the practice 
of computerized tomography, 

2 .1 DMA _COLI£CTjraN ; 

In computerized tomography the number of photons 
transmitted through an object is counted by an array 
of detectors. This array is rotated by a fixed angular 
incremen t between each set of measurements. From these 
measurements the total attenuation along each measurement 
path is derived, which can be written as: 

= exp [ - /|i(r,0) ds ] (2,1) 

= detector reading (counts/second) 

= source strength (counts/second) 

= path of radiation (ray) 

= chord along v\hich s is integrated 
= absorption coefficient 
= cylindrical co-ordinates. 


N 


whe re , 


N 


N 

c 

§ 

c 


r,0 
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From ( 2,1 ) 

jCii(r,0) ds = In (N^M) = P(s,9) (2,2) 


\Ahere, 

P(s,6) is called projection data and it is the 
line integral of function along the line specified by 
s and 0, see Fig, 2.1 , In practical application of 
reconstruction from projections, the measured data 
corresponds to estimates of P(s, ©) for numerous discrete 
values of s and 9, Further, for the simplification of 
the reconstruction algorithm the projections are sampled 
uniformly in both s and 0, Therefore, projections are 
measured at N angles A 9 apart, vdth M equispaced rays 
A s apart for each of the angle of view. 

We choose , 

% 

A © = 

^ ® +ir 

In Older to ensure that the collection of rays 
specified by 

{(mAs, nA0) ; “M^m< + M, 1 < n < N} 
cover the v\hole object from a complete range of directions. 
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Source 



FIG. 2.1 PARALLEL -BEAM GEOMETRY. 



6 


We denote by P(R, © ) the 1-D Fourier transform 
with respect to s of projection at ar^le ©^, where, 0^ 
denotes the nA0 

P(R, 0^) = / PCs, 0^) ds (2.3) 

— .00 


2.2 I%GE RECONSTRUCTION 
PROJECT IONT ' 


BY CONVOLUTION AND BACK- 


The starting point for the derivation of this 
method is the projection theorem for Fourier transform 
or sometimes called projection-slice theorem [6], vtiich 
gives the result 


P(R, 0) =1i(R Cos©, R Sin©) (2.4) 


The physical interpretation of this result is as 
follows s Consider © a fixed angle, then the theorem says 
that as R varies, the value of P(R, 0^) is the same as 

7 V j n 

the value of at R and arrgle © in Fourier space. 

Taking the inverse Fourier transform of (2,4) 

^(x,y) = 5 7 "PCR.®) Cose^ 

O -00 

(2.5) 

For the computation of the integral (involved in 
eqn,(2,5) , v/e need to make some approximation because 
the limits on R varies from -oo to +oo '/\tiich introduces 
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divergence. Therefore, an appropriate cut-off of R is 
used to evaluate the integral, which otherwise diverges. 
So, a function W(R) is introduced in eqn,(2,5) called 
as window function or filter. Window function is chosen 
in such a way that 


W(R) 


Some function of R, [R] ^ R. 

0 , I R I > R^ 


v\iiere, 


R = 1/2 As, As is distance between rays, 

V/indow function of the above form causes abrupt 
truncation of the function in the Fourier domain, which 
leads to oscillations in the spatial domain. The amplitude 
of these oscillations may be reduced greatly by choosing 
a different window function whose value tapers less abruptly 
to zero at the cutoff frequency R^, 


To reduce the amplitude of oscillations in spatial 
doiiain, which arises because of abrupt truncation in 
Fourier domain we developed new family of filters here at 
II T Ka ipur , 


The base filters used in this thesis are? 

Lewitt's filter W(R) =1 ~e|RI/R , iRj < R^ 

c — c 

^0 , iRj > R^ 
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Sin(7rR/R ) 
c 


Sine filter , W(R) = . -ruR/R 


^ R 


R > R. 


Sine square filter , IV(R) = 


Sin(7iR/R.) 2 

■] . R < R 

■C' ^ 




0 


, R > R. 


After introdueing the window function in eqn,(2,5) 


we get 


ti(x,y) = / /^?(R,G) .W(R) ''' {RldR d< 


,0 


( 2 . 6 ) 

Substituting for P(R,0) from eqn,(2,3) and inter- 
changing the order of integration over s and R we get 


It 00 

(x,y) — / / P(s,0) q(x Cos9 + y Sin0 - s) ds d0 

O -00 

(2.7) 

If we assume object to be of unit radius then 
limits on s will be from -1 to +1 , so we get 

% 1 

u.(x,y) = I I P(s,0) q(x Cos© + y Sin0 - s) ds d© 

-o -1 

( 2 . 8 ) 

where. 
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q(s) 


+1/2AS 

I 

-1/2 As 


IRI V'(R) 


_i2nfts 

e 


dR 


(2.9) 


q(s) is called convolving function. To get the 
corresponding higher order filters [7]" of the filter 
functions stitted ^bove we use tho rxjiation, 


q 


( s) 







( 2 . 10 ) 


The equations (2,8) and (2,9) form the basis of 
the method of convolution-back projection. To see how 
this name ( convolution-backprojecti on) arises, we break 
eqn, (2,8) into two equations; 

P(s’,0) = J P(s,©) q(s’-s) ds (2.11) 

r1 

ti(5<,y) = / P(x Cos© y Sin9,©)d© (2.12) 

o 


The operation represented by (2,1 l) is convolution 
and P(s’, ©) called convolved projection at angle © is 
the result of the convolution over s of the projection 
at angle 0 and the convolving fu ration q(s). 

The operation represented by (2,12) is knov/n as 
backpro jecti on and geometrical representation of ihis 
is that the argument of ^ in eqn, (2,10) are parameters 
of -aray through the point (x,y) at ar^le ©, so that p,(x,y) 
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is obtained by the integration of the values of the 
convolved projections associated with all rays passing 
throu^ the point (x,y). 

For the implementation of eqns,(2,11) and (2,12) 
we discretise these equations. The simplest way is to 
use trapezoidal rule and we get. 




iM CO 

2 P(s’, ©„) 

n=1 


(2.13) 


P(m*As, ©j^) ~ AS 


+M 

2 

m=-M 


P(mAs,©j^) q((m‘-m)As) 

(2,14) 


Up till now we v\ere discussing, about parallel-beam 
geometry and this algorithm can also be used for fan-beam 
geometry. In fan-beam geometry each ray is considered 
as one of a set of diverging rays ('OjP), v^iere p determines 
the source position and a is the angle of divergence of the 
ray f rom the source -to -center lir« and line integral or 
projection data is denoted by g(a,p). We assume that source 
is always outside the object to be reconstructed so that 
g(ar,p) is zero for 1 ol greater than some acute angle §, 

From Fig, (2,2), we get the following relationships: 

s = D Sin a 
© = p+ a 

r Cos - (©-0) — s = U Sin( a* -a) 


(2.15a) 

(2.15b) 

(2,15c) 
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= [ r Cos(p- 0) ]^ + [D + r Sii!i(p- 02 ^ 


a 


t 


tan 


[- 


^ Cosif ~ $) ^ 

D + r Siiatp- 0; 


] 


(2.15d) 

(2,15e) 


where, 

D is the redius of the circle on which source lies 
U is the distance between any point ( r, 0) and 
the source 

is the angle of the diverging ray v\hich passes 
through the point (r, 0), 

Writing equivalent of eqn, (2«5) in cylindrical 

co-ordinates 

. 271 ; 00 

p(r,0) = ^ S S P(R,©) exp[i 2idl r Cos(0 - 0)] 
o -00 

IRI dR d0 (2.16) 


Incorporating the window function in eqn, (2,16) and 
substituting P(R, 9) from eqn, (2,3) 


(r,0) 



27t 

I P(s,0) W(R) exp[i2rcR [r Cos(9 - 0)-s]] 
o 

de ds dR 


(2,17) 

Nov; changing variables from ( s, 9) to (a,p) using 
transformations (2,15a,b,c,d,e) we obtain; 
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FIG. 2.2 DIRECT INVERSE GEOMETRY . 
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-1-R 

G 

+6 

2% 

/ 

/ 

X 

-Rc 

r-6 

0 


W(R) ]Rl dp dcrdR (2,18) 

In simple form we can write eqn, (2,18) as; 


D CO 

U (r,0) ^ ^ (cr‘ ,p ) dp 


(2,19) 


v\iiere. 


gc(a',p) 


+ 6 

/ 9( cr, p) Cos(cr) q(U Sin( a’ -a))dcr 

-6 


(2,20) 


and 


+R 


q(s) = / IR|W(R) exp (i 2-n:Rs) dR (2,21-^ 


Discretised version of eqns, (2,19) and (2,20) using 
trapezoidal rule can be written as : 


ti(r,0) = 


D ceAp 
“' 2 ” 


where 


*g ( k a , mAp) 
c 


N-1 . 

2 g (cr*, mAp) 

n=0 L^(p . r,0) 

^ ( 2 , 22 ) 


+M 

S Cos(ia) g(ia» mZip) q(Sin( k*-i)a) 
i=-M 


(2.23) 


So, for the implementation of algorithm we use 
eqns, (2,13) and (2,14) for parallel-beam gecmetry and 
eqns, (2,22) and (2,23) for fan-bean geometry. 
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2.3 DISPLAY; 

In general image is displaved in the form of a 
matrix. The elenents of the matrix vary from 0 to 99 
or from 0 to 33 and are called ‘CT numbers’ , At the 
bottom of the picture are indicated the maximum and 
minimum values of the parameter ‘LITF’ , 'LITF’ is the 
reconstructed value and is related to point densi ty values 
or it is the density of the unit cell if we divide- -the 
entire output matrix into nxn square grid (n can be 
21 or 41), Since we aie interested in cr oss-secti onally 
averaged density. Thus <LITF> is represented as cross 
sectional average density, and it is the sum of values 
of LITF over the entire cross-section of the object 
divided by the number of unit cells present (See Figs. 2, 4, 
2.5). 

We have also represented the reconstructed functions 
by the hidden line 3-diroon§ional grapii using tho graphics 
facility of DEC iogo (See Fig, 2.3), 
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FIG. 2. 3 


HIDDEN LINE 3-DIMENSIONAL REPRESEN- 
TATION OF RECONSTRUCTED FUNCTION 
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CH^TER 

The present study involved the following; 

(a) Implementation and testing of the existing algorithms 
and their application in the point density measurement 
of air~water flow, and 

(b) Implementation of the convolution— backprojection 
algorithm [t] with higher order filters and its 
eventual validation against simulated data for various 
assumed distributions. 


3,1 


Msyjj S JAfl TH. Ejg STJNG FI LTERS; 


The algorithm for fan-beam geometry written in 
FORTRAN was tested for simulated data sets for different 
known distributions. Reconstruction of the known distri- 
butions was done for different value of epsilon for the 
filter suggested by Lewitt [2] . The known distributions 


taken are; 
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n (r) 

= 6 (0) 
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The siinulated projection data for the above cases 
is generated by taking object radius to be unity and 
ray spacing 1,25*^, There are 41 rays in one scan and 
41 different source positions. The reconstructed density 
distribution for above cases is shown in Figs, 3,1, 3,2a, 
3,2b, 3,3, 3,4. The reconstructed values for all the 
cases are in good agreement, within + 0,03 of the true 
values. It is seen that as the value of the epsilon 
goes frail 0 to 1 the value of the average of the distri- 
bution function decreases by 0,01 to 0,03, This is 
expected because of the nature of V/(R), The results are 
summarised in Table 5 in Appendix B, 

The air-water flow data was processed and calibrated 
to get the density distribution for various cases of<P>, 
the cross-sectionally averaged density. Figure 3,5 shows 
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the results for two extreme cases of epsilon values. In 
Fig, 3,6, we have compared the <Pqj , the reconstructed 
density with <P> for various sets (scans)|th^ reconstructed 

3 

density is within + 0,03 g/cm of the actual density 
values. These results are comparable with the results of 
two other methods [8, 9], 


3.2 


JES-TiM 


The algorithm for parallel beam geometry was tested 
for the same known distributions, vhich were used for fan- 
beam geometry. VVhile generating the data ray spacing is 
taken as 0.1 and there are 21 rays in one scan and 21 
source positions. The data generated is summarised in 
Table 1 and Table 2 of Appendix A, 

Then, different filters and their corresponding 
higher order filters were taken and functions were recons- 
tructed. The results are shown in the Fig, 3,11 to Fig, 3,48. 
for the known distribution functions mentioned in the 
previous page. It is seen that the reconstructed values are 
within + 0,01 to + 0,025 of the actual values for most of 
the cases. 


It is seen in all the cases that higher order filters 
give more accurate results and degree of smoothness also 
increases. The results for all the filters and their correspon 
ing higher order filters are sumraa3rf.sed in Table 1 -4 in 
Appendix B, 
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Fig. 3.1 RESULTS WITH SIMULATED DATA. 
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Fig. 3.2(b) Results with simulated data for Lewitt's 
filter with different epsilon values. 
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Fig. 3.3 Results with simulated data for p{r) = 8(0) 
for Lewitts filter with different epsilon 
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Fig. 3.4 Results with simulated data. 
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FIG 3.8 THE SPATIAL DOMAIN REPRESENTATION OF THE 
CONVOLVING FUNCTION WITH LEWITT'S FILTER <6-1.0) 
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FIG 3.9 THE SPATIAL DOMAIN REPRESENTATION OF THE 
CONVOLVING FUNCTION WITH SINC FILTER 
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FIG 3.11 RESULTS WITH SIMULATED DATA F0R;i(r)-1.0 

FOR LEWITT'S I. II. AND III ORDER FILTERS (£«0.0) 
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FIG 3. 12 RESULTS WITH SIMULATED DATA FOR JU(r) 
FOR SING I. II. AND III ORDER FILTERS 
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FIG 3.13 RESULTS WITH SIMULATED DATA FOR ;Ll<r> -1. 0 
FOR SINCSQUARE I. 1 1. AND III ORDER FILTERS 
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FIG 3.14 RESULTS WITH SIMULATED DATA FOR JU<r) -8 (0) 
FOR LEWITT'S I. II. AND III ORDER FILTERS WITH 6-0.0 
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FIG 3.15 RESULTS WITH SIMULATED DATA FOR JLl(r) -S (0) 
FOR SINC I. 1 1. AND III ORDER FILTERS 
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FIG 3.16 RESULTS WITH SIMULATED DATA FOR ^UCr) -6(0) 
FOR SINCSQUARE I, II. AND III ORDER FILTERS 
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FIG 3.18 RESULTS WITH SIMULATED DATA FDR LEWITT'S 
I. II. AND III ORDER FILTERS WITH E -0.0 
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FIG3-/9 RESULTS WITH SIMULATED DATA FDR 
SING I. II. AND III DRDER FILTERS 
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FIG 3.20 RESULTS WITH SIMULATED DATA FOR 
SINC I, II. AND III ORDER FILTERS 
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FIG 3.21 RESULTS WITH SIMULATED DATA FOR 

SINCSQUARE I, II. AND III ORDER FILTERS 
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FIG 3.23 RESULTS WITH SIMULATED DATA FOR 

FOR LEWITT'S I. II. AND III ORDER FILTERS WITH 6-0.0 
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FIG 3.24 RESULTS WITH SIMULATED DATA FOR 
FOR SING I. II. AND III ORDER FILTERS 
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FIG 3.25 RESULTS WITH SIMULATED DATA FDR JUl<r) -r 

FDR SINCSQUARE I. II. AND III DRDER FILTERS 
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FIG 3.26 RESULTS WITH SIMULATED DATA FOR M<f'> “exp <''> 
FOR LEWITT'S I, II. AND III ORDER FILTERS WITH G-0. 0 
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FIG 3.27 RESULTS WITH SIMULATED DATA FOR JUL(r) -exp <r) 
FOR SINC I. 1 1. AND III ORDER FILTERS 
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FIG 3.28 RESULTS WITH SIMULATED DATA FOR JiCr) -axp (r) 
FOR SINCSQUARE I. II. AND III ORDER FILTERS 
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FIG 3.33 RESULTS WITH SIMULATED DATA FOR H<r>- 
FOR LEWITT'S I. II, AND III ORDER FILTERS WITH E 
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FIG 3.34 RESULTS WITH SIMULATED DATA FORJU(r)-r 
FOR LEWITT'S I. II. AND III ORDER FILTERS WITH 6- 
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FIG 3.35 RESULTS WITH SIMULATED DATA FOR 
FOR LEWITT'S I. II, AND III ORDER FILTERS WITH 6 
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FIG 3.36 RESULTS WITH SIMULATED DATA FOR - 
FOR LEWITT'S I. II. AND III ORDER FILTERS WITH E 



0 . 90 


58 



q: 


x 

a: 

LU 

q: 

Q 

UJ 

Q 

LU 

•— 1 

a 

q: 

Q 


tt: 

Q 

q: 

3 

o 


o 

Q CD 


Q 


LU i-i 


2 

Q 

2: on 

in 

a 

cc 

3 h- 

q: 

u 


if) if) 

1— 1 

LU 

X 

if) 

Lu 

O) 

1 — 

< a 



0 
« 

1 


P 




H— 

— 

. . . - I - . 

- i- 

f-- 


■ -f - 

f- 


o 

a 

O 

o 

a 

o 

Q 

a 

o 

GO 


(0 

in 



OJ 

•—1 

o 

d 

« 

Q 

« 

o 

« 

o 

a 

• 

o 

Q 

d 

d 



< 

rf 







FIG 3.37 RESULTS WITH SIMULATED DATA F0R;4(r)- 
FOR LEWITT'S I. II. AND III ORDER FILTERS WITH 6 
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FIG 3,38 RESULTS WITH SIMULATED DATA FOR 
LEWITT'S I. II, AND III ORDER FILTERS WITH 8 
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FIG 3.40 RESULTS WITH SIMULATED DATA FOR -exp <r> 
FOR LEWITT'S I. II. AND III ORDER FILTERS WITH S-l.O 
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FIG 3. 41 RESULTS WITH SIMULATED DATA FOR J-lCr*) “Qxp 
FOR LEWITT'S I. II, AND III ORDER FILTERS WITH 6-1.0 
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FIG 3.42 RESULTS WITH SIMULATED DATA FDR j.l(r) ■£ (D> 
FOR LEWITT'S I. II. AND III ORDER FILTERS WITH 6-1.0 
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FIG. 3.43 HIDDEN LINE 3-DIMENSIONAL REPRESENT- 
ATION OF RECONSTRUCTED FUNCTION 
p(r) =1.0 WITH LEWITT'S FILTER (6=0.0) 



65 



■|G. 3,44 HIDDEN LINE 3-DIMENSIONAL REPRESEN- 
TATION OF RECONSTRUCTED FUNCTION 
Jj(r)=6(0) WITH LEWITT'S FILTER (6 = 0.0) 
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5.3.45 HIDDEN LINE 3-DIMENSIONAL REPRESEN- 
TATION OF RECONSTRUCTED FUNCTION 
p(r) =exp(tr) WITH LEWITT’S FILTER(6=0.0) 
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FIG. 3.46 HIDDEN LINE 3-DIMENSlONAL REPRESEN - 
TATION OF RECONSTRUCTED FUNCTION 
jj(r) =exp(-r) WITH LEWITT'S FILTER (6=0.0) 
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RG. 3.47 HIDDEN LINE 3-DIMENSIONAL REPRESEN- 
TATION OF RECONSTRUCTED FUNCTION 
p(r)=r WITH LEWITT’S FILTER (6 = 0.0) 
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FIG. 3.48 HIDDEN LINE 3-OIMENSIONAL REPRESEN- 
TATION OF RECONSTRUCTED FUNCTION 
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CHAPTER 4 

The linear filters suggested by Lewdtt [2] appear 
to give a reasonably good reconstii-iction of the test- 
functions and of the radial density profile for the two- 
phase air— water flow data [3,6], The error in the recons- 
truction of the average density, <p>, is approximately, 

+ 0,03 g/cm^, which is comparable with various other 
existing techniquesof two-phase density measurement. 

The higher-order filters show an improvement in the 
performance of the Convolution-bac kprojection algorithm 
for test-functions of various types. The experimental data 
was not processed with the higher-order filters because 
tl^i r relative performance could not be judgad in the 
absence of point-wise density measurement by any alternate 
method. 

The following suggestions are expected to improve 
the perf ormance/efficiency of the Convolution Bac kprojection 
algorithm; 

( 1 ) Incorporating better integration schemes for the 
evaluation of convolving functions, 

(b) Implementing still hi^er order filters, i,e. taking 
with k 2, 

(c) Investigating filters other than the linear filters 
tised in the present work. 
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